1. Загрузите датасет life_expectancy_data.RDS (лежит в папке домашнего задания). Это данные с основными показателями, через которые высчитывается ожидаемая продолжительности жизни по метрике World Development Indicator на уровне стран. В данных оставлены строки, относящиеся к положению женщин в 2019 г.
data <- read_rds("C:/Users/79214/Documents/BioStat_2023/GitHub/life_expectancy_data.RDS")
summary(data)
##    Country               Year         Gender          Life expectancy
##  Length:195         Min.   :2019   Length:195         Min.   :55.49  
##  Class :character   1st Qu.:2019   Class :character   1st Qu.:70.02  
##  Mode  :character   Median :2019   Mode  :character   Median :77.55  
##                     Mean   :2019                      Mean   :75.52  
##                     3rd Qu.:2019                      3rd Qu.:80.95  
##                     Max.   :2019                      Max.   :88.10  
##   Unemployment    Infant Mortality      GDP                 GNI           
##  Min.   : 0.178   Min.   : 1.40    Min.   :1.884e+08   Min.   :3.754e+08  
##  1st Qu.: 3.735   1st Qu.: 5.35    1st Qu.:1.117e+10   1st Qu.:1.094e+10  
##  Median : 5.960   Median :13.50    Median :3.967e+10   Median :4.009e+10  
##  Mean   : 8.597   Mean   :19.61    Mean   :4.660e+11   Mean   :4.864e+11  
##  3rd Qu.:10.958   3rd Qu.:30.23    3rd Qu.:2.476e+11   3rd Qu.:2.457e+11  
##  Max.   :36.442   Max.   :75.80    Max.   :2.143e+13   Max.   :2.171e+13  
##  Clean fuels and cooking technologies   Per Capita      
##  Min.   :  0.00                       Min.   :   228.2  
##  1st Qu.: 34.50                       1st Qu.:  2165.3  
##  Median : 80.70                       Median :  6624.8  
##  Mean   : 65.98                       Mean   : 16821.0  
##  3rd Qu.:100.00                       3rd Qu.: 19439.7  
##  Max.   :100.00                       Max.   :175813.9  
##  Mortality caused by road traffic injury Tuberculosis Incidence
##  Min.   : 0.00                           Min.   :  0.0         
##  1st Qu.: 8.20                           1st Qu.: 12.0         
##  Median :16.00                           Median : 46.0         
##  Mean   :17.06                           Mean   :103.8         
##  3rd Qu.:24.00                           3rd Qu.:138.5         
##  Max.   :64.60                           Max.   :654.0         
##  DPT Immunization HepB3 Immunization Measles Immunization Hospital beds   
##  Min.   :35.00    Min.   :35.00      Min.   :37.00        Min.   : 0.200  
##  1st Qu.:85.69    1st Qu.:81.31      1st Qu.:84.85        1st Qu.: 1.301  
##  Median :92.00    Median :91.00      Median :92.00        Median : 2.570  
##  Mean   :87.99    Mean   :86.76      Mean   :87.31        Mean   : 2.997  
##  3rd Qu.:97.00    3rd Qu.:96.00      3rd Qu.:96.50        3rd Qu.: 3.773  
##  Max.   :99.00    Max.   :99.00      Max.   :99.00        Max.   :13.710  
##  Basic sanitation services Tuberculosis treatment Urban population
##  Min.   :  8.632           Min.   :  0.00         Min.   : 13.25  
##  1st Qu.: 62.919           1st Qu.: 73.00         1st Qu.: 41.92  
##  Median : 91.144           Median : 82.00         Median : 58.76  
##  Mean   : 77.380           Mean   : 77.57         Mean   : 59.12  
##  3rd Qu.: 98.582           3rd Qu.: 88.00         3rd Qu.: 78.02  
##  Max.   :100.000           Max.   :100.00         Max.   :100.00  
##  Rural population Non-communicable Mortality  Sucide Rate        continent 
##  Min.   : 0.00    Min.   : 4.40              Min.   : 0.300   Africa  :52  
##  1st Qu.:21.98    1st Qu.:11.85              1st Qu.: 2.050   Americas:38  
##  Median :41.24    Median :17.20              Median : 3.500   Asia    :42  
##  Mean   :40.88    Mean   :17.05              Mean   : 4.802   Europe  :48  
##  3rd Qu.:58.08    3rd Qu.:22.10              3rd Qu.: 6.600   Oceania :15  
##  Max.   :86.75    Max.   :43.70              Max.   :30.100
  1. Сделайте интерактивный plotly график любых двух нумерических колонок. Раскрасть по колонке континента, на котором расположена страна
plot_ly(
  data = data[(data$Unemployment != 0) & (data$`Tuberculosis Incidence` != 0)],
  x = ~ Unemployment,
  y = ~ `Tuberculosis Incidence`,
  color = ~ continent) %>%
    layout(
    title = 'Отношение уровня заболеваемости туберкулезом к уровню безработицы',
    yaxis = list(title = 'Количество заболеваний туберкулезом',
                 zeroline = FALSE),  # Уберём выделения нулевых осей по y
    xaxis = list(title = 'Уровень безработицы',
                 zeroline = FALSE)) # Уберём выделения нулевых осей по y
  1. Проведите тест, на сравнение распределений колонки Life expectancy между группами стран Африки и Америки. Вид статистического теста определите самостоятельно. Визуализируйте результат через библиотеку rstatix.
filter_data <- data %>%
  filter(.$continent == 'Africa' | .$continent == 'Americas') %>%
  select(Country, `Life expectancy`, continent)

filter_data %>%
  ggqqplot(x = 'Life expectancy', facet.by = "continent")

#сравним распределение ожидаемой продолжительности жизни в Америке и Африке, для этого применим критерий Манна-Уитни
filter_data_aremicas <- filter_data %>%
  filter(.$continent == 'Americas')
filter_data_africa <- filter_data %>%
  filter(.$continent == 'Africa')

stat.test <-  filter_data %>%
  wilcox_test(`Life expectancy` ~ continent) %>%
  add_xy_position(x = "continent")
  
stat.test #p-value << 0.05, значит можно отвергнуть нулевую гипотезу, следовательно ожидаемая продолжительность жизни разная в группах Африки и Америки
## # A tibble: 1 × 11
##   .y.       group1 group2    n1    n2 statistic        p y.position groups  xmin
##   <chr>     <chr>  <chr>  <int> <int>     <dbl>    <dbl>      <dbl> <name> <dbl>
## 1 Life exp… Africa Ameri…    52    38       107 6.34e-13       86.6 <chr>      1
## # ℹ 1 more variable: xmax <dbl>
ggboxplot(
  filter_data, 
  x = "continent", y = 'Life expectancy' , 
  ylab = 'Life expectancy', xlab = "Continent", 
  add = "jitter"
  ) + 
  labs(subtitle = get_test_label(stat.test, detailed = TRUE)) + 
  stat_pvalue_manual(stat.test, tip.length = 0)

  1. Сделайте новый датафрейм, в котором оставите все численные колонки кроме Year. Сделайте корреляционный анализ этих данных. Постройте два любых типа графиков для визуализации корреляций.
new_data <- data %>%
  select(-c(Country, Year, Gender, continent))

cor_new_data <- cor(new_data)
corrplot(cor_new_data, method = "number", type = "lower", tl.cex = 0.5, number.cex = 0.5)

ggpairs(new_data,
        title = 'Correlations in Life expectancy dataset',progress = F) +
    theme_minimal() +
    scale_fill_manual(values = c('#69b3a2')) +
    scale_colour_manual(values = c('#69b3a2'))

  1. Постройте иерархическую кластеризацию на этом датафрейме.
scaled_new_data <- scale(new_data)

scaled_clear_dist <- dist(scaled_new_data, 
                      method = "euclidean")

as.matrix(scaled_clear_dist)[1:6,1:6]
##          1        2        3        4        5        6
## 1 0.000000 7.605708 6.331840 4.414874 6.645623 7.923487
## 2 7.605708 0.000000 2.624659 7.921597 3.357361 3.631018
## 3 6.331840 2.624659 0.000000 6.321666 4.350331 3.464837
## 4 4.414874 7.921597 6.321666 0.000000 8.095849 7.161240
## 5 6.645623 3.357361 4.350331 8.095849 0.000000 4.966244
## 6 7.923487 3.631018 3.464837 7.161240 4.966244 0.000000
data_clear_hc <- hclust(d = scaled_clear_dist, 
                        method = "ward.D2")
fviz_dend(data_clear_hc, 
          cex = 0.1)
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

  1. Сделайте одновременный график heatmap и иерархической кластеризации. Содержательно интерпретируйте результат.
pheatmap(scaled_new_data, 
         show_rownames = FALSE, 
         clustering_distance_rows = scaled_clear_dist,
         clustering_method = "ward.D2", 
         cutree_rows = 5,
         cutree_cols = length(colnames(scaled_new_data)),
         angle_col = 45, 
         main = "Dendrograms for clustering rows and columns with heatmap")

#кластеризация помогает разбить колонки примерно на 5 групп, сами описательные признаки образуют 4 группы (т.е 4 региона), причем сильная корреляция по столбцам GDP и GNI в одном регионе. Между собой связны столбцы о безработице, количестве заболеваний туберкулезом, детская смертность, в друггую группы выделены столбцы 'Immunization', отдельная группа - количество мест в больницах и количество суицидов. По первой группе есть корреляция по одному региону.
  1. Проведите PCA анализ на этих данных. Проинтерпретируйте результат.
data_full.pca <- prcomp(new_data, 
                        scale = T)
summary(data_full.pca)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6    PC7
## Standard deviation     2.7526 1.4841 1.3952 1.17177 1.08375 0.96347 0.9288
## Proportion of Variance 0.3988 0.1159 0.1025 0.07227 0.06182 0.04886 0.0454
## Cumulative Proportion  0.3988 0.5147 0.6172 0.68945 0.75126 0.80012 0.8455
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.85740 0.69263 0.68937 0.59106 0.54986 0.47085 0.36596
## Proportion of Variance 0.03869 0.02525 0.02501 0.01839 0.01591 0.01167 0.00705
## Cumulative Proportion  0.88421 0.90946 0.93447 0.95286 0.96877 0.98044 0.98749
##                           PC15    PC16    PC17    PC18      PC19
## Standard deviation     0.34546 0.26941 0.20224 0.06968 1.017e-15
## Proportion of Variance 0.00628 0.00382 0.00215 0.00026 0.000e+00
## Cumulative Proportion  0.99377 0.99759 0.99974 1.00000 1.000e+00
fviz_eig(data_full.pca, addlabels = T, ylim = c(0, 40))

fviz_pca_var(data_full.pca, col.var = "contrib")

fviz_pca_var(data_full.pca, 
             select.var = list(contrib = 3), 
             col.var = "contrib") #очень сильная корреляция переменной 'Life expectations' с 1 компонентой

fviz_contrib(data_full.pca, choice = "var", axes = 1, top = 24)

fviz_contrib(data_full.pca, choice = "var", axes = 2, top = 24)

fviz_contrib(data_full.pca, choice = "var", axes = 3, top = 24)

#первые две компоненты объясняют 50% дисперсии компонент
#большой вклад в анализируемые данные компоненты вносят переменные: Unempoyment, sucide rate, tuberculosis treatment, hospital beds и т.д
#выделяются группы: "Immunization", 1 четверть и граница 2 и 3 четверти.
  1. Постройте biplot график для PCA. Раскрасьте его по значениям континентов. Переведите его в plotly. Желательно, чтобы при наведении на точку, вы могли видеть название страны.
ggbiplot(data_full.pca, 
         scale=0, alpha = 0.1) + 
  theme_minimal()

data_clear_with_ch <- data %>% 
  select(-c(Country, Year, Gender))

ggplotly(ggbiplot(data_full.pca, 
         scale=0, 
         groups = as.factor(data_clear_with_ch$continent), 
         ellipse = T,
         alpha = 0.2) +
  theme_minimal())
  1. Дайте содержательную интерпретацию PCA анализу.
#дополнения к предыдущим объяснениям. Данные могу быть объяснены 3 переменными. Есть переменные, которые имеют отрицательную корреляцию. На последнем графике видны выбросы, данные можно кластеризовать по континентам, но это происходит не очень эффективно.
  1. Сравните результаты отображения точек между алгоритмами PCA и UMAP.
umap_prep <- recipe(~., data = new_data) %>% 
  step_normalize(all_predictors()) %>% 
  step_umap(all_predictors()) %>%  
  prep() %>%  
  juice() 

umap_prep %>%
  ggplot(aes(UMAP1, UMAP2)) + 
  geom_point(aes(color = as.character(data_clear_with_ch$continent)), 
             alpha = 0.7, size = 2) +
  labs(color = NULL)

#есть сходства в результатах отображения точек: точки Африканского континета аггрегируются в левом верхнем углу, как и в PCA, видна аггрегация Европейского континента (хотя оно не такое "сконцентрированное" как в PCA), нет таких сильных выбросов в UMAP. В UMAP точки находятся более плотно, лучше видны кластеры.
  1. Давайте самостоятельно увидим, что снижение размерности – это группа методов, славящаяся своей неустойчивостью. Удалите 5 случайных колонок. Проведите PCA анализ. Повторите результат 3 раза. Наблюдаете ли вы изменения в куммулятивном проценте объяснённой вариации? В итоговом представлении данных на биплотах? С чем связаны изменения между тремя PCA?
#1
rand_1 <- sample(1:19, 5)
new_data_1 <- subset(new_data, select = -c(rand_1))

data_full.pca_1 <- prcomp(new_data_1, 
                        scale = T)
summary(data_full.pca_1)
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.3501 1.3777 1.15047 1.05237 0.97405 0.88640 0.79875
## Proportion of Variance 0.3945 0.1356 0.09454 0.07911 0.06777 0.05612 0.04557
## Cumulative Proportion  0.3945 0.5301 0.62460 0.70370 0.77147 0.82759 0.87317
##                           PC8     PC9    PC10    PC11   PC12    PC13    PC14
## Standard deviation     0.6919 0.63884 0.55333 0.54901 0.3608 0.33158 0.20260
## Proportion of Variance 0.0342 0.02915 0.02187 0.02153 0.0093 0.00785 0.00293
## Cumulative Proportion  0.9074 0.93652 0.95838 0.97991 0.9892 0.99707 1.00000
fviz_eig(data_full.pca_1, addlabels = T, ylim = c(0, 40))

fviz_pca_var(data_full.pca_1, col.var = "contrib")

fviz_pca_var(data_full.pca_1, 
             select.var = list(contrib = 3), 
             col.var = "contrib") #объяснено 55,5% объясненной вариации

graph_1 <- ggplotly(ggbiplot(data_full.pca_1, 
         scale=0, 
         groups = as.factor(data_clear_with_ch$continent), 
         ellipse = T,
         alpha = 0.2) +
  theme_minimal())

#2
rand_2 <- sample(1:19, 5)
new_data_2 <- subset(new_data, select = -c(rand_2))

data_full.pca_2 <- prcomp(new_data_2, 
                        scale = T)
summary(data_full.pca_2)
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.4795 1.2080 1.14397 1.01342 0.96305 0.90564 0.75882
## Proportion of Variance 0.4391 0.1042 0.09348 0.07336 0.06625 0.05858 0.04113
## Cumulative Proportion  0.4391 0.5434 0.63683 0.71019 0.77644 0.83502 0.87615
##                            PC8     PC9    PC10    PC11   PC12    PC13    PC14
## Standard deviation     0.68138 0.61773 0.58059 0.47287 0.3871 0.32157 0.27205
## Proportion of Variance 0.03316 0.02726 0.02408 0.01597 0.0107 0.00739 0.00529
## Cumulative Proportion  0.90932 0.93657 0.96065 0.97662 0.9873 0.99471 1.00000
fviz_eig(data_full.pca_2, addlabels = T, ylim = c(0, 40))

fviz_pca_var(data_full.pca_2, col.var = "contrib")

fviz_pca_var(data_full.pca_2, 
             select.var = list(contrib = 3), 
             col.var = "contrib") #объяснено 54,1% объясненной вариации

graph_2 <- ggplotly(ggbiplot(data_full.pca_2, 
         scale=0, 
         groups = as.factor(data_clear_with_ch$continent), 
         ellipse = T,
         alpha = 0.2) +
  theme_minimal())

#3
rand_3 <- sample(1:19, 5)
new_data_3 <- subset(new_data, select = -c(rand_3))

data_full.pca_3 <- prcomp(new_data_3, 
                        scale = T)
summary(data_full.pca_3)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.3060 1.4008 1.3197 1.12170 1.05576 0.86576 0.76065
## Proportion of Variance 0.3798 0.1402 0.1244 0.08987 0.07962 0.05354 0.04133
## Cumulative Proportion  0.3798 0.5200 0.6444 0.73427 0.81389 0.86743 0.90876
##                            PC8     PC9    PC10   PC11    PC12    PC13      PC14
## Standard deviation     0.69631 0.63187 0.50350 0.3062 0.20290 0.07002 3.191e-16
## Proportion of Variance 0.03463 0.02852 0.01811 0.0067 0.00294 0.00035 0.000e+00
## Cumulative Proportion  0.94339 0.97191 0.99001 0.9967 0.99965 1.00000 1.000e+00
fviz_eig(data_full.pca_3, addlabels = T, ylim = c(0, 40))

fviz_pca_var(data_full.pca_3, col.var = "contrib")

fviz_pca_var(data_full.pca_3, 
             select.var = list(contrib = 3), 
             col.var = "contrib") #объяснено 54,6% объясненной вариации

graph_3 <- ggplotly(ggbiplot(data_full.pca_3, 
         scale=0, 
         groups = as.factor(data_clear_with_ch$continent), 
         ellipse = T,
         alpha = 0.2) +
  theme_minimal())

graph_1
graph_2
graph_3
#при удалении любых 5 столбцов кумулятивный процент объясненной вариации примерно равно 55%, следовательно каждый удаленный столбец поднимает объясненную вариацию примерно на 1%. Эллипсы на биплотах стали более вытяннутыми, т.е данные стали более разбросаны. В зависимости от того, какие столбцы были удалены, расположение эллипсов отличается: эллипсы превращаются в круги, происходят перемещения эллипсов по графику. Все 3 графика достаточно сильно отличаются друг от друга.
  1. Давайте самостоятельно увидим, что снижение размерности – это группа методов, славящаяся своей неустойчивостью. Создайте две дамми-колонки о том: (1) принадлежит ли страна к африканскому континенту, (2) Океании. Проведите PCA вместе с ними, постройте биплоты. Проинтерпрейтируйте результат. Объясните, почему добавление дамми-колонок не совсем корректно в случае PCA нашего типа.
dammi_data <- data %>%
  mutate(Is_African = case_when(.$continent == 'Africa' ~ 1,
                                .$continent != 'Africa' ~ 0),
         Is_Oceanian = case_when(.$continent == 'Oceania' ~ 1,
                                .$continent != 'Oceania' ~ 0))
dammi_data <- dammi_data %>%
  select(-c(Country, Year, Gender, continent))

dammi_data.pca <- prcomp(dammi_data, 
                        scale = T)
summary(dammi_data.pca)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6    PC7
## Standard deviation     2.8400 1.4851 1.3960 1.28878 1.11384 1.01306 0.9514
## Proportion of Variance 0.3841 0.1050 0.0928 0.07909 0.05908 0.04887 0.0431
## Cumulative Proportion  0.3841 0.4891 0.5819 0.66098 0.72006 0.76893 0.8120
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.92891 0.84328 0.70113 0.67720 0.62108 0.55110 0.46766
## Proportion of Variance 0.04109 0.03386 0.02341 0.02184 0.01837 0.01446 0.01041
## Cumulative Proportion  0.85312 0.88699 0.91040 0.93223 0.95060 0.96506 0.97548
##                           PC15    PC16    PC17    PC18    PC19    PC20
## Standard deviation     0.39263 0.35858 0.34148 0.26517 0.20133 0.06900
## Proportion of Variance 0.00734 0.00612 0.00555 0.00335 0.00193 0.00023
## Cumulative Proportion  0.98282 0.98894 0.99449 0.99784 0.99977 1.00000
##                             PC21
## Standard deviation     9.958e-16
## Proportion of Variance 0.000e+00
## Cumulative Proportion  1.000e+00
fviz_eig(dammi_data.pca, addlabels = T, ylim = c(0, 40))

fviz_pca_var(dammi_data.pca, col.var = "contrib")

fviz_pca_var(dammi_data.pca, 
             select.var = list(contrib = 3), 
             col.var = "contrib") #объяснено 55,5% объясненной вариации

ggplotly(ggbiplot(dammi_data.pca, 
         scale=0, 
         groups = as.factor(data_clear_with_ch$continent), 
         ellipse = T,
         alpha = 0.2) +
  theme_minimal())
#добавление двух дамми-колонок уменьшает кумулятивный процент объясненной вариации примерно на 2%, но точки остались на прежних местах, эллипсы остались такими же. То есть наши данные просто стали хуже объясняться глаными компонентами. Это происходит из-за того, что добавленные столбцы бинарные и не несут количественного и качественного смысла без категориальных столцов.